Dynamics of thermalisation in small Hubbard-model systems 
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We study numerically the thermalisation and temporal evolution of a two-site subsystem of a 
fermionic Hubbard model prepared far from equilibrium at a definite energy. Even for very small 
systems near quantum degeneracy, the subsystem can reach a steady state resembling equilibrium. 
This occurs for a non-perturbative coupling between the subsystem and the rest of the lattice where 
relajcation to equilibrium is Gaussian in time, in sharp contrast to perturbative results. We find 
similar results for random couplings, suggesting such behaviour is generic for small systems. 
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Understanding the origin of statistical mechanics from 
a purely quantum-mechanical description is an interest- 
ing area of active research [l|-|9| . Of particular interest is 
the situation of an isolated quantum system partitioned 
into a subsystem and a bath. We ask the question: how 
do observables on the subsystem thermalise when the to- 
tal system is in a pure state? Seminal works [l|, |2| have 
demonstrated the concept of 'canonical typicality' that 
most random pure states of well-defined energy for the 
total system lead to thermalised reduced density matrices 
(RDMs) for the small subsystem. Numerical works have 
demonstrated thermalisation in spin or boson systems for 
various observables of the subsystem [sHsl (and of the en- 
tire closed system [6|,ll0|). Recent theoretical work [11[ 
has investigated whether thermalisation of small subsys- 
tems, initially far from equilibrium, is generic. 

In this letter, we investigate the temporal relaxation 
towards a steady state, focusing on the regime where the 
steady state appears thermalised. We consider a small 
Hubbard ring of fermions away from half filling, with 
two adjacent sites as a subsystem and the other sites as 
a bath. We prepare the system in a product state of sub- 
system and bath pure states in a narrow energy window. 
Even for such a small system, we find a steady-state RDM 
close to a thermal state, down to quantum degenerate 
temperatures. Moreover, we find that the RDM diagonal 
elements approach a steady state as an exponential decay 
for weak subsystem-bath coupling. This becomes a Gaus- 
sian decay at a non-perturbative coupling, with a decay 
rate that departs significantly from the Fermi golden rule. 
We note that this is distinct from the Gaussian behaviour 
in driven systems that remain out of equilibrium J12i ]. 
and in decoherence dynamics [13| of ofF-diagonal RDM 
elements in systems that cannot thermalise. 

The Model. Taking motivation from cold atoms in op- 
tical lattices IJ, ll5| where local addressing is possible, 
we study a local cluster in a generic (non-integrable) 
interacting system with a quasi-continuous spectrum. 
We will examine how a local subsystem {S) thermalises 
with the rest of the system as a bath {B) via unitary 
evolution of the whole system under the Hamiltonian 
H — Hs + Hb + W where Hs {Hb), with eigenstates 



\s)s (|^)b) of energy Es (et), acts solely on the subsystem 
(bath). W is the coupling between the subsystem and 
the bath. For A = 0, the eigenstates are products of sub- 
system and bath eigenstates, denoted \sh), with energies 
Esb — Eg + £()■ The homogeneous case corresponds to 
A = 1. Specifically, we choose a two-site subsystem in an 
L-site Hubbard ring of fermions: 
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where ni„ — c\^Cia is the number operator on site i with 
spin (7. The lattice is a ring with the subsystem sites 
at i = 1,2 and bath sites at z = 3 to L. The hop- 
ping integrals are Ja = J(l + £, sgn{a)), with a non-zero 
^ — 0.05 to remove level degeneracies due to spin rotation 
symmetry. We set the on-site repulsion U = J to give 
us a metallic system with interacting bath modes while 
avoiding the formation of strong features in the many- 
body density of states at U ^ J arising from Hubbard 
interactions. We will let J be the unit of energy. This 
Hamiltonian preserves the total particle number, N, and 
spin component, S^, but not the total spin S^. In this 
work, wc fill the L ~ 9 lattice with 8 fermions of total 
spin S^ — 0. The two-site subsystem has 16 eigenstates 
and the bath has 8281 eigenstates, while the compos- 
ite system has a total of 15876 states with average level 
spacing A ~ 10"'^. This is small enough to allow exact 
diagonalisation, but large enough to provide a smooth 
density of states. 

Consider a system prepared in a pure state of the form 



K 



\^it = 0),Eo)= Y. 

b^=bl 



1 



(2) 



where \si)s is the initial subsystem state, e.g. \ t,t)s 
with parallel spins on the two sites. |^(i = 0)) contains a 
linear combination of B bath eigenstates |fei)_B within an 



energy shell of width 5b, chosen such that (^|if l^t) = Eq. 
The width 6b (= 0.5 in this work) is small on the scale for 
variations in the density of states. The system evolves in 
time: |*(t)) = e~*^*|^(0)). The subsystem is described 
completely by the RDM which traces over the bath states 
|6)b: p{t) = TrB|^(i))(*(i)|. This is evaluated using the 
eigenstates of H from exact diagonalisation. 

Equilibrium States. Before discussing relaxation dy- 
namics during thermalisation, we identify first the pa- 
rameter regime where p does relax to thermal equilib- 
rium. We say that a subsystem thermalises if its RDM 
p{t) approaches the thermal RDM w after some time t 
(shorter than A~^.) The thermal RDM, cj, is diagonal 
with elements Wss — s{s\uj\s)s esc NB{EQ—es,N—ns,S^ — 
sf), where \s)s is a subsystem eigenstate \s)s with en- 
ergy Eg, Us particles and spin sf and Nb{E, rib, sf) is the 
number of bath states with energies in [E, E + SE] with 
rib particles and spin sg . We have to specify energy, num- 
ber and spin because they are globally conserved by the 
Hamiltonian H. We can define an effective temperature 
Tcs = [dlogNB/dE\Eo]~^ provided that the system is in 
a state with energy uncertainty SE ^ A, the level spac- 
ing. (In the thermodynamic limit, uj takes the form of a 
Gibbs canonical distribution [10[ — if particles are not 
exchanged with the bath, utss oc e~'^'''^ .) 
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FIG. 1. (Colour online) Time average of diagonal RDM el- 
ements, Vas, for the ris = 2, Sg = sector as a function 
of composite energy Eq, for initial subsystem states |t4-it)s 
(solid line), |t, t)s (dashed) and \t,i)s (dotted). The four 
elements are labelled by their energies, ascending from ei to 
£4. Thick hues: the corresponding elements for the thermal 
state LJ, found by counting bath states with spin S" — si and 
N — Us particles within a Gaussian energy window of width 
SE — 0.5, centered on energy Eq — Sa- 



We now present our results for a system starting from 
the initial states ([2]) . We avoid the regime of very small 
subsystem-bath coupling A where the subsystem RDM, 
p, is strongly dependent on the initial state even at long 
times due to finite-size effects. Nevertheless, we find that 
even such a small system can reach a steady state for 
couplings A larger than a surprisingly small crossover 



value Ath ^ 1- The RDM becomes virtually diagonal — 
even the sum over the fluctuating off-diagonal elements, 
iJ2s^s' '"ss']^^^' i^ 10^^ to 10^"' smaller than each diago- 
nal element. Fig. [1] shows the steady-state values of the 
diagonal elements of the RDM, r = limi— >.oo /J" p{t)dt/T 
as a function of the composite energy Eq for a coupling 
of A = 0.5. For a variety of initial states, p{t) is only 
weakly dependent on the details of the initial state at 
long times for — 3 < i?o < 6, approaching the thermal 
form cj expected from the canonical ensemble. 




FIG. 2. Dependence on the initial state, Ar (solid), and dis- 
tance from the thermal state, ctui (hollow), as a function of 
coupling A, at composite energies Eo — —2 (circle) and 1.77 
(square). Inset: effective temperature Tch for Eq = —2. 



Next, we establish the range of the coupling A over 
which the system forgets its initial state and thermalises. 
We expect the system to retain memory of the initial 
state at weak coupling (A <C 1). Moreover, for A ^ 1, the 
eigenspectrum becomes significantly altered by the cou- 
pling, splitting into several bands and we see oscillations. 
This is a feature of the projection of the initial state on 
the strongly-coupled link. Therefore, we expect that the 
loss of memory of the initial state and thermalisation are 
possible only in a range of intermediate couplings. To 
quantify this, we calculated the root-mean-square varia- 
tion in diagonal RDM elements due to using different ini- 
tial subsystem states: Ar = ^ X^sK'^ss^) — (rss)^]^ , with 
(. . .) averaging over all 16 initial states in the subsystem 
Fock basis (i.e., eigenstates at J = 0). A small Ar indi- 
cates memory loss. We have also measured the closeness 
to the thermal state uj using a^^ = 5 X]s(kss — ^ss\}- We 
see from Fig. [2] that memory loss and thermalisation oc- 
cur in the intermediate range Ath ^ A < 3 with crossover 
value Ath ^ 0.1 at Eq = -2 and 1.77. 

We also find that the relative probabilities of different 
states in the rig = 2, s^ — sector fit a Boltzmann form: 
logTss — —Es/Tcff + const. For states near the centre of 
the eigenspectrum {E^ ~ 1.77), the effective temperature 
Toff is infinite. At Eq — —2, we find Teff — 2 up to A ~ 2 
(Fig. [2] inset). We estimate the chemical potential to be 
2 J ~ 2 so that, unlike in previous work, we see thermal- 
isation at temperatures down to quantum degeneracy. 



We note that these thermahsed systems are surpris- 
ingly small. Popescu et al. [2] give an estimate of the 
number dji of composite-system eigenstates spanned by 
the initial state sufficient for thermalisation — if the 
probability that CTi^ > y = 0.1 is at least as small as 
X = 0.01, then dn > (97r3/2r2) ln(2/X) ~ 70000. This 
is almost two orders of magnitude larger than the number 
of states (> i5b/A) spanned by our initial state which is 
as low as 950. Moreover, we find thermalisation a.tU — J 
for smaller systems than at t/ <C J. We believe strong 
inelastic scattering in the interacting bath enables effi- 
cient thermalisation a.t U = J when system size is larger 
than the inelastic scattering length (ex ,P /U'^ for small 
U/J). 

Time Evolution. Having established the coupling 
range for thermalisation for model ([l}, we will now dis- 
cuss our main results for the temporal relaxation towards 
the steady state. Fig. [3] shows examples of the time evo- 
lution of the diagonal RDM element pss (t) with s = Si for 
two coupling strengths. The system is again prepared in 
the product state ^ with |sj)s = |tit)s- These results 
are computed for energy Eq = —2. We do not expect 
our results to depend strongly on Eq unless the system 
is close to a strongly correlated ground state. 
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FIG. 4. (Colour online) Rates of decay of ps^si for the subsys- 
tem state Si = It, t) with initial state ^ at Eo = —2 (•/+) 
and 1.77 (B/x) and random model at Eo = -2 (O/A). Weak 
coupling/exponential decay: ~dps-Si/dt\t^o at short times 
found from exponential fits (#,■,0) agree with Fermi golden 
rule prediction, 7fgr (solid lines, gradient 2). Moderate cou- 
pling/Gaussian decay: fit to Gaussian with rate F (x,-l-,A) 
agrees with Fi (dashed lines, gradient 1). Vertical lines mark 
estimates of the crossover values Acxp and Acauss (A^'^ for 
Hubbard/random models). Data in the crossover region are 
rates obtained from attempted fits to either form. 
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in A, the RDM element corresponding to a subsystem 
state s 7^ Si is approximated by p{t): 
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FIG. 3. The RDM element pss as a function of time t with s = 
1 1, t) for the initial state ^ at total energy Eq = -2 (Tcft ~ 
2) with Si = s. Left: coupling A = 0.1 with exponential fit 
(dashed). Right: A = 1 with Gaussian fit (dashed). 



We find qualitatively different relaxation behaviour 
for perturbativc and non-perturbative couplings (Fig. [3]). 
Whereas the RDM relaxes towards the steady state ex- 
ponentially in time at weak coupling (A < Acxp), the 
relaxation follows a Gaussian form at larger coupling 
(A > Acauss)- Interestingly, this Gaussian regime cov- 
ers the coupling range where the system thermalises. 

We can understand our results at short times or 
weak coupling. At short times, we can approximate 
|*(t)) ~ (1 - ii7t)|*(0)). It can be shown (and 
our numerics agree) that the element ps.s.(t) ~ 1 — 
Ffi^ for t < ti = l/max(£'5f, - Eq), with Fi = 
AE^^^^,J(s6|V^|*(0))|2]i/2_ The maximum energy dif- 
ference between states coupled by hopping (V) is of the 
order of the single-particle bandwidth 4 J and so ii ~ 1/4. 

At weak coupling, we can go beyond ti by treating the 
coupling XV as a perturbation to the uncoupled Hamilto- 
nian Hs + Hb- It is readily shown that, to leading order 



after the composite system is prepared in the state ©. 
The element PsiSi is most readily found by using Tr(p) — 
1 to give PsiSi{t) = 1 — J2s^s- Pss{t). This perturbation 
theory is valid until time ^2 when Ps^si has dropped signif- 
icantly below unity. For times between ti and ^2, eq. ([3|) 
follows the Fermi golden rule (FGR): PsiSi{t) decreases 
linearly in time with dps^si/dt oc — A^. Beyond the FGR 
regime, we expect to see exponential decay (see, e.g., ap- 
proximate Markovian schemes of the Lindblad type [16|) 
as is found in our data (Fig.[3]left) for A < Acxp = 0.1. In 
our case, the initial state is not a bath eigenstate. This 
gives small fiuctuations on top of a simple linear-i decay, 
due to interference between terms in the inner sum in ([3]). 

We check in Fig. [4] that the FGR prediction agrees 
quantitatively with dpsiSi/dt\t=o for Eq = —2 and 1.77, 
found from the parameters obtained for the exponential 
fit to ps^s, for t > ti: ps,sdi) ~ Ae-^*"*")/'" + (1 - A). 
The FGR rate, 7fgRj is found by averaging ^ over a 
time t between ti and ^2: —jFGRt = /g [Pstst [t') — l]dt' /t. 
This procedure is needed for a non-zero level spacing A. 
We point out that the exponential fit fails at very weak 
coupling (A ^ 10~^) when the system barely relaxes. 

We now consider larger couplings where A ~ 0(1). In- 



stead of exponential relaxation, we find a good fit (Fig. [3] 
right) to a Gaussian decay: PsiSi{t) ^ Ce~^ * + (1 — C). 
This is seen for couplings A > Acauss = 1- The decay 
rate F now increases linearly with A and is as large as 
the bandwidth scale 1/ii ~ 4. It appears insensitive to 
energy Eq. Interestingly, we see in Fig 4 that, in the 
regime where Ffi '^ 1, the decay rate F is well approxi- 
mated by Fi from the short-time expansion, which sug- 
gests F = Fi/Ci (In our data, 0.97 < C^ < 1.) In other 
words, perturbation theory gives the early-time precursor 
to the full Gaussian form. This suggests the interpreta- 
tion that the time interval of validity of the Fermi golden 
rule (ii < t < 12) narrows and vanishes as A increases to 
Acauss- For coupling range Aoxp < A < Acauss, the be- 
haviour is less clear cut — the decay starts as a Gaussian 
but becomes exponential at later times. The amplitude 
of this exponential tail decreases with increasing A, be- 
coming negligible as A reaches Acauss- 

Random Couplings. To verify that the two regimes of 
relaxation are not specific to our model Hamiltonian, 
we proceeded to study an alternative model where the 
subsystem-bath coupling V is replaced by a random Her- 
mitian matrix W which still respects the conservation of 
the global particle number N and spin S^. Each non-zero 
matrix element of W is Gaussian distributed, with the 
variance chosen such that Tr(VF^) = Tr(y^). Thus, we 
can compare H = Hs+Hb+XV with H = Hs+Hb+\W 
with similar decay rates. In this model, we expect 1/ti 
to be of the order of the full bandwidth ~ 20 for A^ = 8, 
S^ = 0. We find exponential relaxation at weak coupling. 
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0.8, and we recover Gaussian relaxation with 



a linear-A decay rate for A > Aca 



8 (Fig. H hol- 



low symbols). The crossover values, Aexp and Acauss, oc- 
cur at nominally higher couplings than for the Hubbard 
ring ([T|) . They become closer to the Hubbard-ring values 
if we mimic the structure of V by restricting the states 
coupled by W: {s'b'\W\sb) 7^ only if |£;,,b, -^^f,| < 4 J, 
the single-particle bandwidth. 

We summarise our results in Fig. \5\ We have shown 
that a two-site subsystem of the Hubbard model relaxes 
to steady states resembling canonical thermal states, 
even for systems with a handful of sites and at quantum 
degenerate energies. This occurs at a non-perturbative 
coupling between the subsystem and bath, correspond- 
ing to nearly homogeneous systems. In this regime, the 
reduced density matrix p{t) displays Gaussian relaxation 
to the thermal state, with a decay rate F linear in the 
coupling A. This contrasts sharply with the perturba- 
tive regime where p(t) exhibits an exponential relaxation 
with a A^ decay rate. We believe that the Gaussian re- 
laxation to thermalisation is a generic feature of closed 
nanoscale systems, as is supported by our results for ran- 
dom Hamiltonians. 

Finally, we note that it can be shown that Fiii ^ 
XJti ~ A irrespective of system size. The subsystem 
thermalises on the time scale of a few hops between the 



subsystem and the bath, by inelastic collisions of the 
fermions within this timescale. This should be insensitive 
to system size for systems larger than the inelastic scat- 
tering length. Therefore, we speculate that the observed 
Gaussian relaxation should remain for large systems. 
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FIG. 5. Schematic behaviour of the reduced density matrix 
p{t) as a function of subsystem-bath coupling A. Top: relax- 
ation to steady state. Bottom: steady state of p{t). Memory 
of initial state also occurs at large A in the Hubbard case. 
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